set     more off
set     matsize 800
clear   all
global  dpath = "[Userpath]\data and code"

use "$dpath\data\firmdata.dta",clear
g conc=cod/wastewater
sort countycode year

by countycode: egen p70=pctile( conc ) if year<stationyear,p(66.7)
by countycode: egen p70_1=mean(p70) 
drop p70
rename p70_1 p70

by countycode: egen p30=pctile( conc ) if year<stationyear,p(33.3)
by countycode: egen p30_1=mean(p30) 
drop p30
rename p30_1 p30

*High
preserve
drop if conc==.
drop if conc<p70

*1.County-level data

foreach var in output wastewater cod nh3n gas so2 dust {
	bysort countycode year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
	gen l`var' = log(`var')
}
duplicates drop countycode year, force

gen llgdpC = log(lgdpC)
drop if llgdpC == .

gen llgdp = log(lgdp)
drop if llgdp == .

*2. Remove outliers
winsor2 lwastewater, replace cuts(0.5 99.5) trim
winsor2 lcod, replace cuts(0.5 99.5) trim
winsor2 lgas, replace cuts(0.5 99.5) trim
winsor2 lso2, replace cuts(0.5 99.5) trim
winsor2 ldust, replace cuts(0.5 99.5) trim

*3. Generate variables
xtset countycode year

gen manualupstream = 1 if manualyear != 33000
replace manualupstream = 0 if manualupstream == .

gen treated = 1 if autoupstream == 1 & year >= autoyear
replace treated = 0 if treated == .

gen treatedM = 1 if manualupstream == 1 & year >= manualyear
replace treatedM = 0 if treatedM == .

forvalues i=-9(1)-1 {
	local j = -`i'
	gen tN`j'= 1 if year-autoyear == `i' & autoyear != 33000
	replace tN`j' = 0 if tN`j' == .
}

forvalues i=0(1)11 {
	gen t`i' = 1 if year-autoyear == `i' & autoyear != 33000
	replace t`i' = 0 if t`i' == .
}

gen tN5M = tN5+tN6+tN7+tN8+tN9
gen t5M  = t5+t6+t7+t8+t9+t10+t11
gen t9M  = t9+t10+t11
order tN5M,before(tN4) 
order t9M,before(t9) 

cd "$dpath\results"
reghdfe lwastewater tN4-t8 t9M treatedM llgdp, absorb(countycode citycode#year ) 				///
        cluster(citycode)
estimates store Dynamic		
coefplot Dynamic, keep(tN4 tN3 tN2 tN1 t0 t1 t2 t3 t4 t5 t6 t7 t8 t9M) vertical ///
         recast(connect) yline(0) coeflabels(tN4 = "-4" tN3 = "-3" tN2 = "-2"   ///
		 tN1 = "-1" t0 = "0" t1 = "1"  t2 = "2" t3 = "3"  t4 = "4" t5 = "5"     ///
		 t6 = "6"  t7 = "7" t8 = "8"  t9M = "≥9")  								///
		 ytitle("Estimated Coefficients")   									///
		 xtitle("Number of years since the establishment of Automatic Stations") 			///
		 ciopts(recast(rcap) msize(medium))	
graph save "Graph" "Figure IX Upper.gph",replace
restore


*Medium/Low
preserve
drop if conc==.
drop if conc>=p70

*1.County-level data

foreach var in output wastewater cod nh3n gas so2 dust {
	bysort countycode year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
	gen l`var' = log(`var')
}
duplicates drop countycode year, force

gen llgdpC = log(lgdpC)
drop if llgdpC == .

gen llgdp = log(lgdp)
drop if llgdp == .

*2. Remove outliers
winsor2 lwastewater, replace cuts(0.5 99.5) trim
winsor2 lcod, replace cuts(0.5 99.5) trim
winsor2 lgas, replace cuts(0.5 99.5) trim
winsor2 lso2, replace cuts(0.5 99.5) trim
winsor2 ldust, replace cuts(0.5 99.5) trim

*3. Generate variables
xtset countycode year

gen manualupstream = 1 if manualyear != 33000
replace manualupstream = 0 if manualupstream == .

gen treated = 1 if autoupstream == 1 & year >= autoyear
replace treated = 0 if treated == .

gen treatedM = 1 if manualupstream == 1 & year >= manualyear
replace treatedM = 0 if treatedM == .

forvalues i=-9(1)-1 {
	local j = -`i'
	gen tN`j'= 1 if year-autoyear == `i' & autoyear != 33000
	replace tN`j' = 0 if tN`j' == .
}

forvalues i=0(1)11 {
	gen t`i' = 1 if year-autoyear == `i' & autoyear != 33000
	replace t`i' = 0 if t`i' == .
}

gen tN5M = tN5+tN6+tN7+tN8+tN9
gen t5M  = t5+t6+t7+t8+t9+t10+t11
gen t9M  = t9+t10+t11
order tN5M,before(tN4) 
order t9M,before(t9) 

cd "$dpath\results"
reghdfe lwastewater tN4-t8 t9M treatedM llgdp, absorb(countycode citycode#year ) 				///
        cluster(citycode)
estimates store Dynamic		
coefplot Dynamic, keep(tN4 tN3 tN2 tN1 t0 t1 t2 t3 t4 t5 t6 t7 t8 t9M) vertical ///
         recast(connect) yline(0) coeflabels(tN4 = "-4" tN3 = "-3" tN2 = "-2"   ///
		 tN1 = "-1" t0 = "0" t1 = "1"  t2 = "2" t3 = "3"  t4 = "4" t5 = "5"     ///
		 t6 = "6"  t7 = "7" t8 = "8"  t9M = "≥9")  								///
		 ytitle("Estimated Coefficients")   									///
		 xtitle("Number of years since the establishment of Automatic Stations") 			///
		 ciopts(recast(rcap) msize(medium))	
graph save "Graph" "Figure IX Lower.gph",replace
restore